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Context. In previous papers, we have shown that, as the rotation of a neutron star slows down, it will be internally heated as a con- 
sequence of the progressively changing mix of particles (rotochemical heating). In previously studied cases (non-superfluid neutron 
stars or superfluid stars with only modified Urea reactions), this leads to a quasi-steady state in which the star radiates thermal photons 
for a long time, possibly accounting for the ultraviolet radiation observed from the millisecond pulsar J0437-4715. 
Aims. For the first time, we explore the phenomenology of rotochemical heating with direct Urea reactions and uniform and isotropic 
superfluid energy gaps of different sizes. 

Methods. We first do exploratory work by integrating the thermal and chemical evolution equations numerically for different energy 
gaps, which uncovers a rich phenomenology of stable and unstable solutions. To understand these, we perform a stability analysis 
around the quasi-steady state, identifying the characteristic times of growing, decaying, and oscillating solutions. 
Results. For small gaps, the phenomenology is similar to the previously studied cases, in the sense that the solutions quickly converge 
to a quasi-steady state. For large gaps (> 0.05 MeV), these solutions become unstable, leading to a limit-cycle behavior of periodicity 
~ 10 6 ~ 7 yr, in which the star is hot (T s > 10 5 K) for a small fraction of the cycle (~ 5 - 20% ), and cold for a longer time. 

Key words, stars: neutron — dense matter — stars: rotation — pulsars: general — pulsars: individual (PSR J0437-4715) 



1. Introduction 



The observation of thermal emission from the surface of a 
neutron star (NS) has the potential to provide constraints on 
its inner structure. In the existing literature, several detailed 
cooling calculations have been compared to the few estimates 
available for the surface temperatures of neutron stars (see 
lYakovlev & Pefhickll2004l for a review and references). These 
calculations are based on passive cooling, at first neutrino- 
dominated, and later driven by photon emission at ages > 10 5 
yr. 

Several proposed mechanisms could keep NSs hot be- 
yond the standard cooling timescale ~ 10 7 yr (see, e.g., 
ISchaab et alJ {l999), but note that this paper gives an incor- 
rect pararnejeriration_of_rotod^^ 

in IGonzalez & Reiseneggell |2010|) . iGonz alez & Reisenegger 
(1201 Ol) compared the abilities of these mechanisms to reheat 
NSs, concluding that two of them might be important in old NSs, 
name ly, reheating by frictional motion of superflu id neutron vor- 
tices ( Alp ar et al.ll9 84: Shibazaki & Lamb 1989) and rotochem- 
ical heating. The latter was first proposed by Reisene ggeti d 1 9951) 
and then improved by Fernandez & Reisenegger (2005), and 
considered the internal structure of non-superfluid NSs with re- 
alistic equations of state (EOSs) in the framework of general rel- 
ativity. It works as follows. As the NS's rotation rate decreases, 
the reduction in the centrifugal force makes it contract. This per- 
turbs each fluid element, raising the local pressure and causing 
deviations from beta equilibrium. Eventually, the system reaches 
a quasi-steady configuration, where the rate at which spin-down 
modifies the equilibrium concentrations is the same as that at 
which neutrino reactions restore the equilibrium. These reactions 
heat the stellar interior, making the star emit thermal radiation. 



The potential presence of superfluid nucleons in the NS's inte- 
rior has b een widely considered to m odel their thermal evolution 
(see, e.g. lYakovlev & Pe thick 2004 for cooling models) since it 
considerably reduces the neut rino reactions and the specific heat 
involving superfluid species ( Yak ovlev et all 1200 ll) . and opens 
new neutrino emissi on processes, nam ely pair breaking and for- 
mation reactions (Flowers et al., fl976h . 



In a previous paper dPetrovich & Reiseneggeti 1201 Ol) . we 
modeled rotochemical heating of millisecond pulsars with only 
modified Urea reactions in the presence of uniform and isotropic 
Cooper pairing gaps of neutrons A„ and protons A p . We veri - 
fied the order-of-magnitude predictions of Reiseneggei] d!997l) . 
finding that the chemical imbalances in the star grow up to the 
threshold value A f /, r = min(A„ + 3A p , 3A„ + A p ), which is higher 
than in the quasi-steady state achieved in the absence of su- 
perfluidity. Therefore, the old superfluid NSs will take longer 
to reach the quasi-steady state than their nonsuperfluid counter- 
parts, and they have a higher a luminosity in this state, given by 
L ;« s - (1 -4) x 10 32 (A jAr /MeV)(p_ 20 /P ms ) erg s" 1 , where 

P_2o is the period derivative in units of 10~ 2l) and P ms is the 
period in milliseconds. With the previous relation, we found 
that energy gaps in the range 0.05 [MeV] < A thr < 0.45 [MeV] 
are consistent with the ultraviolet emission of PSR J0437-4715 
dKargaltsev et alll2004l) . 

We extend our previous analysis to the case where the much 
faster direct Urea reactions are allowed. We find a qualitatively 
new behavior of rotochemical heating, where the temperature 
and chemical imbalances oscillate around the quasi-steady state. 

The structure of this paper is the following. In Sect. [2] we 
present the basic equations of rotochemical heating and the 



1 Present address: Department of Astrophysical Sciences, 
Princeton University, Princeton, NJ 08544, USA, e-mail: 
cpetrovi@astro . princeton . edu 



2 



Petrovich and Reisenegger: Long-period thermal oscillations in superfluid MSPs. 



phase-space integrals for direct Urea reactions with Cooper pair- 
ing gaps and chemical imbalances. In Sect. [3] we describe our 
results for the thermal evolution and the linear stability analysis 
of old NSs. We summarize our main conclusions in Sect. [4] 



larger during the quasi-steady state, lengthening the timescale to 
arrive at this state compared with the non-superfluid case, and 
predicting higher temperatures in old NSs. In this paper, we in- 
clude the powerful direct Urea reactions to our previous study. 



2. THEORETICAL FRAMEWORK 

2.1. Rotochemical heating: basic equations 

The basi c framework for rotochemical h eating is explained in 
detail in lFernandez & Re isenegger (2005), and the modifications 
made to consider the Cooper p airing effects are described in 
iPetrovich & Reiseneggerl(l2010h . The latter is the framework we 
use throughout this paper. Therefore, we just point out the fun- 
damental equations for completeness and to clarify the notation 
of the present paper. 

We consider the simplest model of a neutron star core, com- 
posed of neutrons, protons, electrons, and muons (npep matter), 
ignoring the potential presence of exotic particles. 

The internal temperature, redshifted to a distant observer, 
Too, is taken to be uniform inside the star because we are model- 
ing the thermal evolution over time scales much longer than the 
diffusion time (Reisene ggerl Il995h . Thus, the evolution of the 
internal temperature for an isothermal interior is given by the 
thermal balance equation (lThorneLll977l) 



rri r CO T cc \ 



(1) 



where C is the total heat capacity of the star, is the total power 
released by the heating mechanism, the total power emitted 
as neutrinos due to Urea reactions, and L™ the power released as 
thermal photons. 

The amount of energy released by each Urca-type reaction is 
rjnpi = Hn -ftp ~Mi ( I - e ' P)> where p-, is the chemical potential of 
the particle species i. Thus, we write the total energy dissipation 
rate as 



r CO OO * 1 1 , CO * 1 i 

L H ~ tlnpe^ 1 npe + Vnpp^ 1 npfi, 



(2) 



where Af np i = f „^ p i - f p /^„ is the net reaction rate of the Urea 
reaction integrated over the core involving the lepton /. 

The photon luminosity is calculated by assuming black- 
body radiation = AncrR^T^ ra , where Rco and r j oo are the 
radius and the surface temperature of the star measured by 
an observer at infinity, respectively. To relate the internal and 
the surface tempera tures, the fully accreted envelope model of 
iPotekhin et all d 19971) is used. 

The evolution of the redshifted chemical imbalances, also 
uniform throughout the core, is given by 



tfnpn ~ Z„pAT npe 



' ZuppATupp + 2.W n pp£l£l, 



(3) 
(4) 



where the terms Z„ p , Z npe , Z np p, 



W npe , and W„ P p are constants 



that depend on the stellar structur e and are kept unchanged with 
respect to their latest definition in Reiseneg ger et alJ d2006h . and 
QQ is the product of the angular velocity and its time derivative 
(propo rtional to the spin-down power). 

As IPetrovich & R eisenegger (2010) showed, the results of 
the evolution with rotochemical heating when computing 
and Ar„ pe in the presence of superfluid nucleons are substan- 
tially different from their s uperfl uid counterparts calculated by 
Fernandez & Reisenegger (2005), since superfluidity strongly 
inhibits these reactions. Thus, the chemical imbalances become 



2.2. Cooper pairing 

In the core, neutrons are believed to form Cooper pairs because 
of their interaction in the triplet 3 P2 states via the anisotropic 
channels \mj\ = (type B) or \mj\ = 2 (type C), while pro- 
tons f orm isotropic, singlet 'So pairs (type A) ( Yak ovlev et all 
2001). Additionally, in the outermost core and inner crust, neu- 
trons are believed to form singlet-state 'So pairs. The 3 P2 (type 
B and C) state description is rather uncertain in the sense that 
the energetically most probable state of nn-pairs Qtnj\ = 0, 1,2) 
is not known, being extremely sensitive to the still unknown nn- 
interaction (see, e.g. Amundsen & 0stgaard 1985). Taking this 
classification into account, Villain & Haensell (2005) solve nu- 
merically the suppression due to each type of superfluidity of 
the net reaction rate for direct Urea and modified Urea reac- 
tions out of beta equilibrium, finding that the suppression due 
to type A superfluidity is of strength between the suppression 
due to anisotropic channels type B and type C superfluidity, re- 
spectively. For simplicity, we consider the energy gaps for the 
neutrons A„ and the protons A ;; at zero temperature, redshifted 
to a distant observer, as parameters that are isotropic ('So pairs) 
and uniform throughout the core of the NS. 

The phase transition for a nucleon species into a superfluid 
state takes place when its temperature falls below a critical 
value T c . This temperature is related to the energy gap at zero 
temperature A(T = 0); for the isotropic pairing channel 'So, 
A(T = 0) = 1.764AT C . Additionally, when the transition occurs, 
the amplitude of the energy gap depends on the temperature by 
means of the BCS equation dYakovlev et al.ll2001l). which can be 
fitted by the practical formula of Levenfish & Yakovlevl (1 19941) 
for the isotropic gap 



6 = 



ACQ 
kT 



= Vi -r/r c |i.456- 



0.157 1.764 
+ 



vttt; t/t c 



(5) 



where 5 is the variable used in the phase-space integrals in Sect. 
2.31 It is straightforward to check that the limiting cases are re- 
produced by Eq. (0, i .e. 6 — when T - T c and 6 = A(T = 
Q)/kT when T « T c .. iLevenfish & Yakovlevl (11994 claim that 
intermediate values of T/T c are also reproduced by this formula 
with a maximum error less than 5%, which is accurate enough 
for the purposes of this work. 

Having defined the energy gap of the nucleon A, with i = 
n, p, it is possible to express the momentum dependence of the 
nucleon en ergy e;(p;) near the Fe rmi level, i.e. [p; - pf.\ <sz pp t , 
as follows ( Yak ovlev et all 1200 ll) 



£i(pi) = pi - ^v 2 F {pi - p F: ) 2 + Aj if pi < p F ., 

€i(pi) = pi + yjv 2 F i P i - p F ,) 2 + A 2 if Pi>p Fi , (6) 

where pp., Vp t , and pi are the momentum, the Fermi momen- 
tum, the Fermi velocity, and the chemical potential of species 
i — n,p, respectively. 

2.3. Neutrino emissivity 

The fastest reactions in NS cores are the direct Urea processes 
n — > p + e~ + v/ (7) 
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p + e 



n + vi, 



where I — e,ji. If these reactions are kinematically allowed, their 
emissivity is m uch greater than those produced by the modified 
Urea reactions ( lYakovlev et alll2001l) . 

We write the neutrino emissivity and the net reaction rate of 
direct Urea reactions involving the lepton / and integrated over 
the core, respectively, as 



AT, 



n pi 



■rlD.TT oc: 



where constants L„i and L p \ are defined in terms of the neutrino 
luminosities for a nonsuperfluid NS in beta equilibrium, as 



jeq r 
Too Jo 



Anr 1 e h S a {n)e^' s 'dr, 



(8) The phase-space integrals Io, e and Io,r cannot be solved 
analyt ically when one o r two n ucle on species are superfluid. 
Thus, [Villain & Haensell d2005l) and iPetrovich & Re isenegger 
d2010h computed these integrals numerically and described 
the suppression produced by the Cooper pairs by means of 
the so-called reductions factors. The latter are defined as 
the ratio of the phase-space integrals with superfluid particle 
species to the non-superfluid integrals, i.e. Rp^ - Id,e/Fd(^i) 
and Rdx - Jz>,r Hereafter, we calculate these inte- 
grals numerically by using the numerical method explained in 
IPetrovich & Re isenegger! ( 20101) . which is based on the Gauss- 
Laguerre quadrature, taking advantage of the exponentially de- 
caying behavior of the Fermi functions appearing in the inte- 
grand. 

3. RESULTS AND DISCUSSION 



(9) 
(10) 



(11) 



where the term 5; is a slowly varying func tion of the baryon 
number density n (e.g.. lYakovlev et ail 200 lh . and A and <t> are 
the usual Schwarzschild metric terms. The quantities Id £ and 
Idx are dimensionless phase-space integrals that contain the de- 
pendence of the emissivity and the net reaction rate, respectively, 
on the chemical imbalances rQ l and the energy gaps A„ and A p . 

To introduce these integrals, it is useful to define the usual 
dimensionless variables normalized by the thermal energy kT, 
as follows 



x = , x v = — , dnci c = , 

j kT kT , «,< kT , 



which represent the energy of the non-superfluid degenerate par- 
ticle j, the neutrino, and the chemical imbalance involving the 
lepton /, respectively, while for the superfluid nucleon i we write 



VF,(Pi-pF,) I , 2 

= — and Zi = sgn(x,) yjxf + 6?, 



where 6, is defined in Eq. ([5]) in terms of A, . Thus, in terms of 
these variables, 



'D,e 



914 
5040tt 6 



and 
Idt = 



914 
5040tt 6 



y^»CO /-»co /-»co /~»CO 

I djCyJCy III dX-yy d X X g 

%JQ %J — CO %J — CO \J — CO 

X f(Zn)f(Zp)f(X e ) [S(x v + £ - Z„ - Zp - X e ) 

+ S(x v - & - z„ - z p - x e )\ (14) 

/"•OO s~*CXl rtCO /"»CO 

I dXyXy ({I ' Xftd X pCl Xq 

%jQ %J —CO \_J — CO \J — CO 

X f(Zn)f(Zp)f(x e ) [d(x v +^l-Zn~Zp- X e ) 



6(x v -%l-Z n -Zp 



4 



where /(■) is the Fermi function f(x) - 1/(1 + e x ), and the nu- 
merical factor in front of the integral is normalizes Id e to 1 when 
the energy gaps and the chemical imbalances are zero. 

In the nonsuperfluid case (i.e. 6„ = 8 D = 0), these integ rals 
reduce to the polynomials calculated by Reiseneggei] (1 19951) 



IdA6„ = 6 P = 0) = F D (&) = 1 + 



1071^, 2 315$ 21$ 



WiSn = S p = 0) = H D (gO 



457tt 2 
714^ 



457tt 4 457tt' 

3 



420^ 42£ 
+ 



457^ 2 457tt 4 457tt 6 



3.1. Evolution 

In this section, we show the peculiar behavior of rotochemical 
heating in MSPs with superfluid nucleons and direct Urea reac- 
tions. To do so, we numerically compute the evolution of Eqs. 
([TJ, (0), and ©. Hereafter, we omit the 00 subscript in the tem- 
perature and the 00 superscript in the chemical imbalances to 
simplify the notation. 

The NS interior stru cture is modeled by the BPAL21 EOS 
of Prakash et al. (1988) fo r the core, supplemented with those 
of P ethick et alJ d 19951) and lHaensel & Pichonl dl994h for the in- 
ner and outer crust, respectively. This model opens electron di- 

(12) re ct Urea reactions in the core when the NS mass is greater than 
1 .67 M and muon direct Urea reactions when its mass is greater 
than 1.69 M , and reaches a maximum mass value of 1.71 M Q . 
Thus, the direct Urea reactions are allowe d within the ma s s rang e 
1.76 + 0.2 M of the MSP J0437-4715 (IVerbiest et all 12008). 
the only MSP whose thermal radiation has been measured so 

(13) far and, therefore , there is an estimate of its surface temperature 
dKargaltsev et alll2004l) . 

For numerical calculations, the evolution of OQ is com- 
puted by assuming magnetic dipole braking with no field decay, 
and relating the magnetic field on the magnetic equator, rota- 
tion period, and period derivative by the conventional formula 
B = 3.2 ■ 10 l9 (Ppyt 2 G, where P is measured in seconds. The 
magnetic field is chosen to be B — 2.8 • 10 8 G, so as to match the 
currently measured values of P and P of PSR J0437-4715. 

In Fig. 1, we show the evolution for three different values 
of the neutron energy gap, with non-superfluid protons and no 
muons, illustrating the appearance of strong thermal oscillations 
at larger gaps. The left top panel, where A„ = 0.01 MeV, shows 
that rotochemical heating converges rapidly to a stable solution 
given by the quasi-steady solution calculated from f = f] npe - 0, 
as it happens when no superflu i d nuc leon species are consid- 
ered dFernandez & Reiseneggerl 120051) or when superfluid nu- 
(15) cleons are present, but only mod ified Urea reactions are allowed 
(IPetrovich & Reisen egger. 2010). A hint of an oscillation is seen 
before the solution converges. 

In the left bottom panel of Fig. 1 , we increase the value of the 
energy gap to A„ = 0.03 MeV. In this case, there is a clearly vis- 
ible damped oscillation before the variables settle to their quasi- 
state values. 

In the upper right panel of Fig. 1, we use A„ = 0.05 MeV 
to show strong oscillations without a clear damping effect. This 
behavior can be understood as follows. Once the chemical imbal- 
ance increases to a critical value ~ A„, the reaction rates increase 
(17) dramatically, suddenly increasing the temperature and thus de- 
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Fig. 1. Evolution of the internal temperature T and the chemical imbalance rj npe for a 1.68 M G star without muons built with the 
BPAL21 EOS (Pra kash et all fl988l) . with the initial conditions T = 10 8 K and n npe = 0. The spin-down is assumed to be caused 
by magnetic dipole radiation, with dipolar field strength B = 2.8 ■ 10 s G and initial period Pq — 1 ms. The dashed lines are the 
quasi-steady state solution calculated from t = f] npe = 0. In all cases, the protons are taken to be non-superfluid, while the neutron 
superfluid energy gaps are A n = 0.01 MeV {upper left panel), A„ = 0.03 MeV (lower left panel), and A„ = 0.05 MeV [upper and 
lower right panels). The lower right panel shows a close-up of the second peak in temperature of the evolution plot shown in the 
upper right panel (A„ = 0.05 MeV) with a linear scale in time. 



creasing the imbalance. At this point, the imbalance is substan- 
tially below the energ y gap and the direct Urea rea ctions are 
strongly inhibited (see Petro vich & ReiseneggerllToiOl for more 
details on the reduction factors). Thus, the star begins to cool 
almost exclusively by photon emission, and the spin-down com- 
pression forces the chemical imbalance to grow almost unim- 
peded by the neutrino reactions until it reaches its critical value 
again. This process is repeated cyclically in a timescale given, 
essentially, by the time taken by the chemical imbalance to come 
back to its critical value ~ A„, where it erases any record of the 
previous cycles. 

The lower right panel of Fig. 1, shows a close-up of one of 
the peaks depicting the evolution with A„ = 0.05 MeV in order 
to resolve and illustrate the rather sharp peaks. We also illustrate 
in this figure that the timescale involved in the rapid increase 
in temperature is on the order of ~ 10 5 yr, while the decrease 
happens on a timescale of ~ 10 6 yr. 

This unexpected behavior will be the focus of our analysis. 

3.2. Stability analysis 

In this section, we analyze the local stability of the temperature 
and chemical imbalance by perturbing the evolution Eqs. ([D and 
(fJJ around the quasi-steady state f = i]„ pe = and computing 
the growth rates of these small displacements. For simplicity, we 
ignore the presence of muons in our analysis. 

We can safely neglect the modified Urea reactions of the 
electrons because the direct Urea reactions are much more ef- 



ficient and the superfluid suppression is similar in both cases. 
Thus, the evolution equations (with j]„ pe in units of temperature) 
can be written as 

T = I [L D [ft/ r - I 6 ] T 6 - L y T a ) = T(T, ij npe ), (18) 
flnpe = ~- {L D I r T 5 ) + ^^QQ = Q{T, r, npe ). (19) 

We considered the fully accreted envelope model of 
Pote khin et all (1 19971) to express the photon luminosity in terms 
of the internal temperature as L y = L y T a , where L y is a numer- 
ical factor that depends on the stellar structure, and a = 2.42. 
Additionally, we defined the functions T and Q to save notation 
in the subsequent analysis. 

We write the perturbed solutions to these equations as 

T = r< s + 5Te yt , and ^ = rj% e + Srj^, (20) 

where T qs and rfi, s pe are the solutions of T — Q — 0, which 
are unique for each combination of input parameters, i.e. en- 
ergy gaps, spin-down power, and structure constants. We define 
y as the growth rate of the perturbations ST and 6j]„ pe . Replacing 
these perturbed solutions in the differential equations leads to the 
usual eigenvalue problem for the growth rates of these perturba- 
tions 

( ^- j) UlH (21) 
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where I is the 2x2 identity matrix and J is the Jacobian ma- 
trix resulting from the linearization of Eqs. ( TT8l > and ( TT9l > and 
evaluated at T qs and n q lpe . Thus, we write the Jacobian matrix el- 
ements, replacing the quasi-steady state and taking logarithmic 
derivatives, as 



J11 

Jl2 
J21 
J22 



ar = Ly Ta _, 



dT 



c 

Ly 



d\nT 



+ 6 



+ (6 - a) 



dT] npe C \ d\nr/ npe 

dg = 2w„ pe nn ( din(zj r ) 

dT ~ 



kT \ d\nT 
2W npe nn /<91n(£J r ) 



dr], 



npe 



hi 



npe 



d In 



T]npe 



(22) 
(23) 
(24) 
(25) 
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Fig. 2. Real and imaginary parts of the growth rates y + and 
y_ (upper panel) calculated from Eq. d26l i and components of 
the Jacobian matrix (lower panel) in Eqs. d22l> . d23l , d24l i. and 
d25l ). for different values of A„. In both panels, the stellar pa- 
rameters are the same as in Fig. 1 and we assume that QQ = 
-3 • l(T 9 [rad 2 /s 3 ]. 

We solve for the two growth rates y + and y_ by simply cal- 
culating 



T ± 



r± = 



Vt 2 - 4A 



(26) 



where r and A are the trace and the determinant of J, respec- 
tively. The stability of this system of differential equations, 
at least when the linear approximations is valid, i.e. close to 
quasi-steady solutions, depends on these growth rates. Thus, if 
R e (j±) < 0, the fixed point solutions will be stable, as happens 
when we have only modified Urea reactions or no superfluid nu- 
cleons are considered. Otherwise, if Re(y ± ) > the solutions 
will diverge. Moreover, the solution will oscillate if Im(y ± ) + 0, 
and the oscillations become important during the evolution if 
\Im(y ± )\ are comparable to or greater than \Re(y ± )\. 

In the upper panel of Fig. 2, we show the real and imaginary 
parts of the two growth rates for different values of A„. We fix the 
value of the spin-down parameters to QQ = - 3 ■ 1 0~ 9 [rad 2 / s 3 ] in 
order to compare our results to those of Fig. 1, where QQ reaches 
this value at the age of 5 • 10 6 yr and slowly changes until ~ 10 s 
yr. When the energy gap is A„ = 0.01 MeV (top panel of Fig. 1), 
the system is locally stable because Re(y ± ) < 0, and the oscilla- 
tion is highly damped because \Re(y ± )\ » \Im(-y±)\, Then, if we 
take A„ = 0.03 MeV (middle panel of Fig. 1), we can see that 
the system is still locally stable, but Re(y + ) and Re(yJ) increase 
in such a way that |^e(y±)| < \Im{y ± )\l(2n), and therefore the 
oscillations are less damped than those with lower energy gaps. 
Up to here, the linearized system predicts that the fixed point T qs 
and ril s pe behaves as a stable spiral (also known as an attractor 
or a sink) in the T - rj„ pe space, hence the linear analysis ensures 
that it is also correct f or the non-linear system close to this fixed 
point (Utrogatz, 1993). 

Increasing the energy gap to a value higher than A„ = 0.03 
MeV, the linearized system becomes unstable because Re(y±) 
changes from being negative to positive. Beyond this point, the 
non-linear system acts as an unstable spiral (also known as a 
repeller or source) in the T - r]„ pe space. If the initial condition of 
the system of Eqs. (TT~8T > and ( fT9l l is close to the fixed point, it will 
diverge from it. Afterwards, the amplitude of the perturbations 
becomes large enough to make the non-linear terms significant, 
and the linear analysis is no longer valid. We show in section 
Sect. l3.4l that. for high energy gaps, the highly non-linear system 
converges to a stable limit cycle. 

3.3. Superfluidity driving the oscillations 

The oscillations are exclusively caused by the appearance of the 
superfluid energy gaps at the Fermi level of the nucleons. This 
drastically changes the dependence of the neutrino emissivity 
and the net reaction rate on the internal temperature and chemi- 
cal imbalance, which causes the oscillations, as we show in this 
section. 

In the lower panel of Fig. 2, we show the previous effect by 
means of the components of the Jacobian matrix in Eqs. d22l . 
i23i . i24i . and 1251 . evaluated at the quasi-steady state. From 
here, we can observe that close to the quasi-steady state, the 
evolution represented by Eqs. ( fTST l and ( fT9T > becomes unstable, 
essentially, because the component Jn = dT /dT changes from 
being negative to positive as we increase A„, changing the sign 
of the trace of the Jacobian t. 

To understand why this happens, we analyze the dimen- 
sionless net heating function £Jp - h that appears in Eq. <TT~8T >. 
For the non-superfluid case, one finds from the polynomials 
given in Eqs. ( TToT l and ( TTTl l that, in the quasi-steady state, 
where T « 77, £J r - h ~ 21/(4577r 6 )(?7/r) 6 , which gives 
<91n(£Jr — I e )/d]nT ~ -6 , i.e. from Eq. (T221 J remains neg- 
ative. In contrast, for the superfluid case, we find that this func- 
tion can increase with temperature for values of the energy gaps 
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sufficiently greater than the temperature. We know that both of 
the functions /p and I £ increase with temperature since this in- 
creases the phase-space availab le to the particles. However, as 
|PetrovicrT & Reiseneg gerl d20 1 Ol) showed by means of the reduc- 
tion factors, superfluidity blocks the phase-space of the integrals 
I e and Iy quite differently, an d for the quasi-steady state we fin d 
that R e <«; R r (see Fig. 2 of iPetrovich & Reisen egger (2 010Id . 
which forces g e Ir - h to be an increasing function of the temper- 
ature. Finally, this happens until the temperature is high enough 
to make the reduction factors comparable and the function de- 
creases with temperature in a similar way to its non-superfluid 
counterpart. 

3.4. Limit cycle in j] npe - T: existence and ^-dependence 

In the upper panel of Fig. 3, we show the evolution of rotochem- 
ical heating in the r\ npe - T space for different initial conditions 
that converge to the same the limit cycle. Moreover, we know 
that the trajectories inside the limit cycle are repelled by the 
fixed point of T and Q, which in this case is an unstable spiral, 
as discussed in Sect. 13.21 Thus, we could construct a trapping 
region, excluding the fixed point, such as the vector that field 
points inward at the boundary of this region and, therefore, ar- 
gue by means of the Poincare-Bendixson theo rem that there is a 
closed orbit, i.e. a stable limit cycle (see, e.g. Strogatz 1994|for 
details on limit cycles). 

In the lower panel of Fig. 3, we show different limit cycles 
by changing the energy gap and label three sections of the evo- 
lution curve with A„ = 0.05 MeV, namely A-B, B-C, and C-A. 
Along the path A-B, the chemical imbalance is ~ A„ and some 
reactions are opened, increasing the temperature until it reaches 
node B, without changing the chemical imbalance. Along the 
path B-C, the temperature is high enough to open abruptly many 
reactions that tend to restore the beta equilibrium, and r\ npe drops 
rapidly. Finally, along the path C-A, the chemical imbalance is 
sufficiently below the energy gap to strongly inhibit the neutrino 
reactions. Thus, no heating mechanism is present and the star 
cools by photon emission. Moreover, no restoring mechanism 
for the beta equilibrium is present and the chemical imbalance 
grows by the spin-down compression. 

We observe from this figure that the amplitude of the limit 
cycles increases with the gap A„, which we can easily under- 
stand from the previous analysis. For the path A-B, the chemical 
imbalance must be larger for larger A„ to open reactions. Then, 
from B to C, we note from the linear analysis in Sect. 13.21 that 
for larger gaps the system becomes more unstable or, more pre- 
cisely, the heating function ^ e Ir - I e becomes more sensitive to 
small perturbations in temperature reaching higher temperatures, 
forcing the chemical imbalances to reach lower values in node C. 
Finally, from C to A, all the curves are parallel since the photon 
emission and spin-down compression do not depend on A„. 

3.5. Limit cycle timescales and the MSP J0437 

We can observe from the bottom panel of Fig. 1 that there are 
two very different scales governing the limit cycles. The first and 
shorter one involves the rapid increase in temperature and ends 
with a rapid decrease in the chemical imbalance, i.e. it goes from 
A to C in the lower part of Fig. 3. The second scale goes from 
C to A and is substantially dominant in the rotochemical heating 
evolution. 

The first is quite stable and on the order of the timescales 
shown in Fig. 2 for the growth rates y, i.e. a timescale of ~ 10 5 
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Fig. 3. Evolution of rotochemical heating in the r\ npe - T space 
for a superfluid star with the same parameters of Fig. 1. Upper 
panel: evolution with A„ = 0.05 MeV for different initial con- 
ditions (dotted-dashed lines) that converge to the same limit cy- 
cle (solid line). The quasi-steady state is indicated (star). Lower 
panel: evolution of the limit cycle for A„ = 0.04, 0.05, 0.06 
MeV. The arrows indicate the path of the cycles A-B, B-C, and 
C-A for the star with A„ = 0.05 MeV. 



yr. However, the second scale depends strongly on the energy 
gap, because the latter sets the amplitude of the drop in r] npe . It 
also depends on since it accounts for the rate at which the 
chemical imbalance grows. For the examples displayed in Figs. 
1 and 3, the longer timescale is ~ 10 6 - 10 7 yr. 

Considering the large difference between both timescales, 
we check that the fraction of the time in which the thermal emis- 
sion would be detectable, i.e. the surface temperature T s > 10 5 
K, is in the range ~ 5-25%, mainly depending on the energy 
gaps and spin-down parameters. Moreover, it is also unlikely to 
find it close to its quasi-steady state and therefore are unable to 
draw further conclusion s about the pairing gaps neede d to ex- 
plain observations as in lPetrovich & Reis enegger (i2010l) . 

However, we can account for the maximum energy gap com- 
bination A,/„. = A„ + A p at which rotochemical heating evolves 
to a stable solution, i.e. its quasi-steady state, and predict a sur- 
face temperature. Thus, for the spin-down paramete rs of MSP 
J0437 and the set of EOSs from lPrakash et al.1 d!988l) that allow 
direct Urea rea ctions of electrons a nd muons for its mass range 
1.76 ± 0.2 M dDelleret all 120081) . namely BPAL21, BPAL31, 
BPAL 32, and BPAL33, the maximum gap combination that pre- 
dict a quasi-steady state is in the range of A,/ lr ~ 0.01 -0.1 MeV. 
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Thus, the temperatures predicted when including superfluidity 
are a few times higher than those predic ted by the non-superfluid 
case dFernandez & Reiseneggerl 12005). They are quite close to 
10 5 K and for some cases can explain the lik ely thermal emission 
of the MSP J0437 dKargaltsev et all l200l . 

Finally, given the rapid increase in temperature, we could ask 
whether the finite thermal diffusion timescales are relevant to our 
study or make any difference to the validity of our assumption of 
an isothermal NS interior. We claim that this assumption is fairly 
safe, since the increase in temperature takes ~ 10 5 yr, while the 
thermal diffusion time is just a few decades for these low tem- 
peratures. 




log t [yr] 




log t [yr] 

Fig. 4. Evolution of the internal temperature T and the chemi- 
cal imbalances r\ npe and t] npjl with the initial condition T — 10 8 
K and null chemical imbalances. The spin-down is assumed to 
be caused by the magnetic dipole radiation, with dipolar field 
strength B = 2.8 ■ 10 8 G and initial period of Pq — 1 ms. 
The dashed lines indicate the energy gap A„ = 0.05 MeV. 
Uppe r panel: evolut i on of a 1.69 M star built with the BPAL21 
EOS (Pra kash et all 1 1988b . which allows muon direct Urea re- 
actions, and we add the vortex creep reheating term of Eq. d27l > 
with an excess of angular momentum J - 10 43 erg s. Lower 
panel: evolution of a 1.68 M star built with the BPAL21 EOS 
dPrakash et all Il988l) . which forbids direct Urea reactions with 
muons. 



3.6. Effect of muons and other reheating mechanisms. 

The presence of muons does not introduce changes to the pre- 
vious analysis, as long as direct Urea reactions with muons are 
allowed. On the other hand, if they react only by means of the 
modified Urea channels, the oscillations are damped and finally 
vanish. 

In the upper panel of Fig. 4, we observe that the only ef- 
fect to be considered is that the oscillations would be driven 
by the reactions involving muons, instead of those with elec- 
trons. This happens because the departure from beta equilibrium 
with muons, i.e. rj npfl , is more strongly affected by the spin-down 
compression, and the direct Urea reactions with muons will lead 
the oscillations since they reach the critical value ~ A„ first. 

Moreover, including another reheating mechanism does not 
alter the existence of the oscillations as long as the temper- 
ature increase produced by this mechanism is not so high as 
to overcome the quasi- steady value of the temperature. As 
iGonzalez & Reiseneggerl (120101) argued, the two most relevant 
mechanisms to reheat a n old neutron star are rotochemical 
heatin g and vortex creep (lAlpar et al.lll984l IShibazaki & L amb 
Il989l) . Therefore we include the latter in the upper panel of Fig. 
4 to show that the oscillations persist even when another reheat- 
ing mechanism is included. The energy dissipation r ate that must 
be inc luded as a heating term in Eq. ([T} is given by dAlpar et all 
1984) 



= 70, 



(27) 



where J is the excess of angular momentum. By choosing 
an intermediate value of J = 10 43 erg s for the models in 
Shibazaki & Lamb ( 1989), we observe that the only effect of the 
extra heating term is to provide a higher floor to the internal tem- 
perature. This results in a straightforward way from the analysis 
in Sect. 13.21 since the Jacobian element in Eq. d22l i is unaltered 
by the new term in Eq. d27l l. However, if this mechanism keeps 
the NS at substantially higher temperatures, it could displace the 
quasi-steady state and change the behavior of the oscillations. 

In the lower panel of Fig. 4, we show the evolution when di- 
rect Urea reactions with muons are forbidden, but modified Urea 
reactions involving muons are present. At first, the chemical im- 
balance of the muons reaches a value close to the energy gaps 
and reheats the star by means of modified Urea reactions. The 
chemical imbalance of the electrons then reaches this threshold 
value and the oscillations start to work normally until the chemi- 
cal imbalance of the muons reaches its quasi-steady state, which 
predicts higher temperatures than those without muons (see the 
bottom panel of Fig. 1). Hence, the oscillations vanish and the 
quasi-steady state is mostly governed by the chemical imbalance 
of the muons. 



3. 7. Effect on classical pulsars 

In this section, we show the effects of these oscillations on clas- 
sical pulsars. In Fig. 5, we illustrate their unstable behavior by 
showing the evolution of rotochemical heating for a NS with an 
energy gap A„ = 0.05 MeV, electron direct Urea reactions, no 
muons, and two magnetic fields of B — 10 11 G and B = 10 12 G. 

From this figure, we observe that the temperature and the 
chemical imbalance can oscillate only during half of the period 
for B = 10 11 G or a full period for B — 10 12 G since the rotation 
rate is insufficient for the imbalance to recover from its drop and, 
therefore, produce any other oscillation. Thus, after the sudden 
increase in temperature at ~ 10 5 yr, the star cools almost exclu- 
sively by means of photon emission, while the chemical imbal- 
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Fig. 5. Evolution of the internal temperature T and the chemical 
imbalance rj npe for a 1.68 M Q star built with the EOS BPAL21 
and no muons, with the initial condition T = 10 9 K and null 
chemical imbalance. The spin-down is assumed to be caused by 
the magnetic dipole radiation, with dipolar field strengths B = 
10 11 G (dashed line) and B = 10 12 G (solid line) initial period 
of P () = 5 ms. The dotted-dashed line indicates the energy gap 
A„ = 0.05 MeV. 

ance grows slowly, never reaching the threshold of the energy 
gap again. 

If we changed the initial period to lower values (e.g. Po = 1 
ms), we could observe from the case with B — 10 11 G one full 
oscillation period and another half, similar to the lower panel of 
Fig. 4. Similarly, lower values of the energy gap could produce 
more oscillations. In contrast, if we increase the initial period, 
the rotation energy to increase the chemical imbalance may not 
be enough to reach the energy gap threshold, and oscillation may 
not occur at all. 

In conclusion, the oscillating behavior generally does not 
persist in classical pulsars, and the standard cooling governs 
their thermal evolution. 
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4. Conclusions 

We have studied a new kind of oscillations of the temperature 
and chemical imbalances, produced during the evolution of ro- 
tochemical heating when the direct Urea reactions are allowed 
in superfluid MSPs with relatively high uniform and isotropic 
Cooper pairing gaps. 

The oscillations work as follows. The direct Urea reactions 
are strongly blocked until the chemical imbalances grow up to a 
threshold value A,/,, = A„ + A p . Soon after this happens, strong 
reactions are suddenly turned on, increasing the temperature and 
decreasing the chemical imbalances. As this happens, the reac- 
tions are again strongly blocked and the star cools down by pho- 
ton emission, while the chemical imbalances again increase until 
they reach the threshold value required to repeat the cycle after 
10 6-7 yr. The temperature stays high for only a small fraction of 
the cycle, making it difficult to detect this effect and predict a 
NS temperatures at a given time. For gaps below ~ 0.05 MeV 
or with muons reacting only by modified Urea processes, the os- 
cillations vanish and the system reaches a quasi-steady state, as 
previously found in the non-superfluid or modified Urea cases. 
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